<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
 <head>
  <title>slagts.f</title>
 <meta name="generator" content="emacs 21.3.1; htmlfontify 0.20">
<style type="text/css"><!-- 
body { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.default   { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.default a { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
span.string   { color: rgb(188, 143, 143);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.string a { color: rgb(188, 143, 143);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
span.comment   { color: rgb(178, 34, 34);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.comment a { color: rgb(178, 34, 34);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
 --></style>

 </head>
  <body>

<pre>
      SUBROUTINE <a name="SLAGTS.1"></a><a href="slagts.f.html#SLAGTS.1">SLAGTS</a>( JOB, N, A, B, C, D, IN, Y, TOL, INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  -- LAPACK auxiliary routine (version 3.1) --
</span><span class="comment">*</span><span class="comment">     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
</span><span class="comment">*</span><span class="comment">     November 2006
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Scalar Arguments ..
</span>      INTEGER            INFO, JOB, N
      REAL               TOL
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Array Arguments ..
</span>      INTEGER            IN( * )
      REAL               A( * ), B( * ), C( * ), D( * ), Y( * )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Purpose
</span><span class="comment">*</span><span class="comment">  =======
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  <a name="SLAGTS.19"></a><a href="slagts.f.html#SLAGTS.1">SLAGTS</a> may be used to solve one of the systems of equations
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     (T - lambda*I)*x = y   or   (T - lambda*I)'*x = y,
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  where T is an n by n tridiagonal matrix, for x, following the
</span><span class="comment">*</span><span class="comment">  factorization of (T - lambda*I) as
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     (T - lambda*I) = P*L*U ,
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  by routine <a name="SLAGTF.28"></a><a href="slagtf.f.html#SLAGTF.1">SLAGTF</a>. The choice of equation to be solved is
</span><span class="comment">*</span><span class="comment">  controlled by the argument JOB, and in each case there is an option
</span><span class="comment">*</span><span class="comment">  to perturb zero or very small diagonal elements of U, this option
</span><span class="comment">*</span><span class="comment">  being intended for use in applications such as inverse iteration.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Arguments
</span><span class="comment">*</span><span class="comment">  =========
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  JOB     (input) INTEGER
</span><span class="comment">*</span><span class="comment">          Specifies the job to be performed by <a name="SLAGTS.37"></a><a href="slagts.f.html#SLAGTS.1">SLAGTS</a> as follows:
</span><span class="comment">*</span><span class="comment">          =  1: The equations  (T - lambda*I)x = y  are to be solved,
</span><span class="comment">*</span><span class="comment">                but diagonal elements of U are not to be perturbed.
</span><span class="comment">*</span><span class="comment">          = -1: The equations  (T - lambda*I)x = y  are to be solved
</span><span class="comment">*</span><span class="comment">                and, if overflow would otherwise occur, the diagonal
</span><span class="comment">*</span><span class="comment">                elements of U are to be perturbed. See argument TOL
</span><span class="comment">*</span><span class="comment">                below.
</span><span class="comment">*</span><span class="comment">          =  2: The equations  (T - lambda*I)'x = y  are to be solved,
</span><span class="comment">*</span><span class="comment">                but diagonal elements of U are not to be perturbed.
</span><span class="comment">*</span><span class="comment">          = -2: The equations  (T - lambda*I)'x = y  are to be solved
</span><span class="comment">*</span><span class="comment">                and, if overflow would otherwise occur, the diagonal
</span><span class="comment">*</span><span class="comment">                elements of U are to be perturbed. See argument TOL
</span><span class="comment">*</span><span class="comment">                below.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  N       (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The order of the matrix T.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  A       (input) REAL array, dimension (N)
</span><span class="comment">*</span><span class="comment">          On entry, A must contain the diagonal elements of U as
</span><span class="comment">*</span><span class="comment">          returned from <a name="SLAGTF.56"></a><a href="slagtf.f.html#SLAGTF.1">SLAGTF</a>.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  B       (input) REAL array, dimension (N-1)
</span><span class="comment">*</span><span class="comment">          On entry, B must contain the first super-diagonal elements of
</span><span class="comment">*</span><span class="comment">          U as returned from <a name="SLAGTF.60"></a><a href="slagtf.f.html#SLAGTF.1">SLAGTF</a>.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  C       (input) REAL array, dimension (N-1)
</span><span class="comment">*</span><span class="comment">          On entry, C must contain the sub-diagonal elements of L as
</span><span class="comment">*</span><span class="comment">          returned from <a name="SLAGTF.64"></a><a href="slagtf.f.html#SLAGTF.1">SLAGTF</a>.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  D       (input) REAL array, dimension (N-2)
</span><span class="comment">*</span><span class="comment">          On entry, D must contain the second super-diagonal elements
</span><span class="comment">*</span><span class="comment">          of U as returned from <a name="SLAGTF.68"></a><a href="slagtf.f.html#SLAGTF.1">SLAGTF</a>.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  IN      (input) INTEGER array, dimension (N)
</span><span class="comment">*</span><span class="comment">          On entry, IN must contain details of the matrix P as returned
</span><span class="comment">*</span><span class="comment">          from <a name="SLAGTF.72"></a><a href="slagtf.f.html#SLAGTF.1">SLAGTF</a>.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Y       (input/output) REAL array, dimension (N)
</span><span class="comment">*</span><span class="comment">          On entry, the right hand side vector y.
</span><span class="comment">*</span><span class="comment">          On exit, Y is overwritten by the solution vector x.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  TOL     (input/output) REAL
</span><span class="comment">*</span><span class="comment">          On entry, with  JOB .lt. 0, TOL should be the minimum
</span><span class="comment">*</span><span class="comment">          perturbation to be made to very small diagonal elements of U.
</span><span class="comment">*</span><span class="comment">          TOL should normally be chosen as about eps*norm(U), where eps
</span><span class="comment">*</span><span class="comment">          is the relative machine precision, but if TOL is supplied as
</span><span class="comment">*</span><span class="comment">          non-positive, then it is reset to eps*max( abs( u(i,j) ) ).
</span><span class="comment">*</span><span class="comment">          If  JOB .gt. 0  then TOL is not referenced.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">          On exit, TOL is changed as described above, only if TOL is
</span><span class="comment">*</span><span class="comment">          non-positive on entry. Otherwise TOL is unchanged.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  INFO    (output) INTEGER
</span><span class="comment">*</span><span class="comment">          = 0   : successful exit
</span><span class="comment">*</span><span class="comment">          .lt. 0: if INFO = -i, the i-th argument had an illegal value
</span><span class="comment">*</span><span class="comment">          .gt. 0: overflow would occur when computing the INFO(th)
</span><span class="comment">*</span><span class="comment">                  element of the solution vector x. This can only occur
</span><span class="comment">*</span><span class="comment">                  when JOB is supplied as positive and either means
</span><span class="comment">*</span><span class="comment">                  that a diagonal element of U is very small, or that
</span><span class="comment">*</span><span class="comment">                  the elements of the right-hand side vector y are very
</span><span class="comment">*</span><span class="comment">                  large.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  =====================================================================
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Parameters ..
</span>      REAL               ONE, ZERO
      PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      INTEGER            K
      REAL               ABSAK, AK, BIGNUM, EPS, PERT, SFMIN, TEMP
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          ABS, MAX, SIGN
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      REAL               <a name="SLAMCH.113"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>
      EXTERNAL           <a name="SLAMCH.114"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           <a name="XERBLA.117"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Executable Statements ..
</span><span class="comment">*</span><span class="comment">
</span>      INFO = 0
      IF( ( ABS( JOB ).GT.2 ) .OR. ( JOB.EQ.0 ) ) THEN
         INFO = -1
      ELSE IF( N.LT.0 ) THEN
         INFO = -2
      END IF
      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.128"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="SLAGTS.128"></a><a href="slagts.f.html#SLAGTS.1">SLAGTS</a>'</span>, -INFO )
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( N.EQ.0 )
     $   RETURN
<span class="comment">*</span><span class="comment">
</span>      EPS = <a name="SLAMCH.135"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'Epsilon'</span> )
      SFMIN = <a name="SLAMCH.136"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'Safe minimum'</span> )
      BIGNUM = ONE / SFMIN
<span class="comment">*</span><span class="comment">
</span>      IF( JOB.LT.0 ) THEN
         IF( TOL.LE.ZERO ) THEN
            TOL = ABS( A( 1 ) )
            IF( N.GT.1 )
     $         TOL = MAX( TOL, ABS( A( 2 ) ), ABS( B( 1 ) ) )
            DO 10 K = 3, N
               TOL = MAX( TOL, ABS( A( K ) ), ABS( B( K-1 ) ),
     $               ABS( D( K-2 ) ) )
   10       CONTINUE
            TOL = TOL*EPS
            IF( TOL.EQ.ZERO )
     $         TOL = EPS
         END IF
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( ABS( JOB ).EQ.1 ) THEN
         DO 20 K = 2, N
            IF( IN( K-1 ).EQ.0 ) THEN
               Y( K ) = Y( K ) - C( K-1 )*Y( K-1 )
            ELSE
               TEMP = Y( K-1 )
               Y( K-1 ) = Y( K )
               Y( K ) = TEMP - C( K-1 )*Y( K )
            END IF
   20    CONTINUE
         IF( JOB.EQ.1 ) THEN
            DO 30 K = N, 1, -1
               IF( K.LE.N-2 ) THEN
                  TEMP = Y( K ) - B( K )*Y( K+1 ) - D( K )*Y( K+2 )
               ELSE IF( K.EQ.N-1 ) THEN
                  TEMP = Y( K ) - B( K )*Y( K+1 )
               ELSE
                  TEMP = Y( K )
               END IF
               AK = A( K )
               ABSAK = ABS( AK )
               IF( ABSAK.LT.ONE ) THEN
                  IF( ABSAK.LT.SFMIN ) THEN
                     IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK )
     $                    THEN
                        INFO = K
                        RETURN
                     ELSE
                        TEMP = TEMP*BIGNUM
                        AK = AK*BIGNUM
                     END IF
                  ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN
                     INFO = K
                     RETURN
                  END IF
               END IF
               Y( K ) = TEMP / AK
   30       CONTINUE
         ELSE
            DO 50 K = N, 1, -1
               IF( K.LE.N-2 ) THEN
                  TEMP = Y( K ) - B( K )*Y( K+1 ) - D( K )*Y( K+2 )
               ELSE IF( K.EQ.N-1 ) THEN
                  TEMP = Y( K ) - B( K )*Y( K+1 )
               ELSE
                  TEMP = Y( K )
               END IF
               AK = A( K )
               PERT = SIGN( TOL, AK )
   40          CONTINUE
               ABSAK = ABS( AK )
               IF( ABSAK.LT.ONE ) THEN
                  IF( ABSAK.LT.SFMIN ) THEN
                     IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK )
     $                    THEN
                        AK = AK + PERT
                        PERT = 2*PERT
                        GO TO 40
                     ELSE
                        TEMP = TEMP*BIGNUM
                        AK = AK*BIGNUM
                     END IF
                  ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN
                     AK = AK + PERT
                     PERT = 2*PERT
                     GO TO 40
                  END IF
               END IF
               Y( K ) = TEMP / AK
   50       CONTINUE
         END IF
      ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Come to here if  JOB = 2 or -2
</span><span class="comment">*</span><span class="comment">
</span>         IF( JOB.EQ.2 ) THEN
            DO 60 K = 1, N
               IF( K.GE.3 ) THEN
                  TEMP = Y( K ) - B( K-1 )*Y( K-1 ) - D( K-2 )*Y( K-2 )
               ELSE IF( K.EQ.2 ) THEN
                  TEMP = Y( K ) - B( K-1 )*Y( K-1 )
               ELSE
                  TEMP = Y( K )
               END IF
               AK = A( K )
               ABSAK = ABS( AK )
               IF( ABSAK.LT.ONE ) THEN
                  IF( ABSAK.LT.SFMIN ) THEN
                     IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK )
     $                    THEN
                        INFO = K
                        RETURN
                     ELSE
                        TEMP = TEMP*BIGNUM
                        AK = AK*BIGNUM
                     END IF
                  ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN
                     INFO = K
                     RETURN
                  END IF
               END IF
               Y( K ) = TEMP / AK
   60       CONTINUE
         ELSE
            DO 80 K = 1, N
               IF( K.GE.3 ) THEN
                  TEMP = Y( K ) - B( K-1 )*Y( K-1 ) - D( K-2 )*Y( K-2 )
               ELSE IF( K.EQ.2 ) THEN
                  TEMP = Y( K ) - B( K-1 )*Y( K-1 )
               ELSE
                  TEMP = Y( K )
               END IF
               AK = A( K )
               PERT = SIGN( TOL, AK )
   70          CONTINUE
               ABSAK = ABS( AK )
               IF( ABSAK.LT.ONE ) THEN
                  IF( ABSAK.LT.SFMIN ) THEN
                     IF( ABSAK.EQ.ZERO .OR. ABS( TEMP )*SFMIN.GT.ABSAK )
     $                    THEN
                        AK = AK + PERT
                        PERT = 2*PERT
                        GO TO 70
                     ELSE
                        TEMP = TEMP*BIGNUM
                        AK = AK*BIGNUM
                     END IF
                  ELSE IF( ABS( TEMP ).GT.ABSAK*BIGNUM ) THEN
                     AK = AK + PERT
                     PERT = 2*PERT
                     GO TO 70
                  END IF
               END IF
               Y( K ) = TEMP / AK
   80       CONTINUE
         END IF
<span class="comment">*</span><span class="comment">
</span>         DO 90 K = N, 2, -1
            IF( IN( K-1 ).EQ.0 ) THEN
               Y( K-1 ) = Y( K-1 ) - C( K-1 )*Y( K )
            ELSE
               TEMP = Y( K-1 )
               Y( K-1 ) = Y( K )
               Y( K ) = TEMP - C( K-1 )*Y( K )
            END IF
   90    CONTINUE
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="SLAGTS.302"></a><a href="slagts.f.html#SLAGTS.1">SLAGTS</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>
